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A  Complete  Mode  Sum  for 
LF,  VLF,  ELF  Terrestrial  Radio  Wave  Fields' 

J.  Ralph  Johler  and  Leslie  A.  Berry 

Contribution  From  the  Central  Radio  Propagation  Laboratory, 
National  Bureau  of  Standards,  Boulder,  Colo. 

Previous  papers  about  VLF,  ELF  propagation  in  the  earth-ionosphere  waveguide  have  emphasized 
the  characteristics  of  the  individual  modes.  In  this  paper,  emphasis  is  on  the  total  field  found  by  sum- 
ming all  modes  vt^hich  contribute  appreciably.  Solutions  are  given  and  calculations  are  made  for  a 
homogeneous  isotropic  ionosphere,  a  homogeneous  ionosphere  with  a  constant  radial  magnetic  field 
imposed,  and  an  ionosphere  with  a  magnetic  field  of  arbitrary  dip  angle. 

In  the  isotropic  case,  one  mode  is  sufficient  to  describe  the  field  at  ELF,  but  at  the  upper  end  of 
the  VLF  band,  two  or  more  modes  may  propagate  7500  km.  When  a  magnetic  field  is  imposed  on  the 
ionosphere,  three  abnormal  modes  or  fields  are  generated  at  the  anisotropic  boundary.  These  abnor- 
mal modes  have  little  effect  on  the  total  field  in  the  VLF  band,  but  they  change  the  ELF  field  appre- 
ciably. Indeed,  in  the  absorption  region  between  VLF  and  ELF,  the  abnormal  modes  may  become 
dominant  beyond  1000  km. 

1.  Introduction 

The  object  of  this  analysis  is  a  more  accurate  description  of  the  field  of  the  LF,  VLF,  ELF  radio 
wave  in  the  space  between  the  ionosphere  and  the  terrestrial  sphere.  A  number  of  investigators 
have  devised  theoretical  methods  to  this  end;  indeed,  considerable  emphasis  has  been  placed  on 
the  development  of  field-strength  curves,  both  as  a  function  of  distance  and  frequency.  Wait  and 
Murphy  [1957]  pubhshed  field-strength  curves  computed  with  the  aid  of  the  geometric-optical 
theory,  applicable  to  short  distances,  and  Wait  [1957,  1962]  and  Wait  and  Carter  [1960]  developed 
curves  with  the  aid  of  the  residue  methods,  or  "mode  theory,"  which  is  a  full-wave  solution  for  an 
assumed  model,  but  employed  the  Debye  [1910]  approximation  for  the  spherical-wave  function. 
This  approximation  was  used  previously  by  March  [1912]  and  von  Rybczynski  [1913]  for  the 
groundwave  problem  and  is,  as  pointed  out  by  March  [1912],  not  valid  in  areas  corresponding  to  the 
most  important  terms  in  the  residue  series  at  frequencies  greater  than  8  or  10  kc/s.  It  is  interest- 
ing to  note  that  von  Rybczynski  [1913]  analyzed  and  corrected  March's  [1912]  paper  in  the  use  of 
this  approximation.  Budden  [1961a]  considered  the  fuU-wave  solutions  or  "mode  characteristics" 
assuming  a  flat  model  for  the  terrestrial  sphere  and  ionosphere,  with  the  imphcation  that  this 
procedure  simphfied  the  analysis.  Alpert  [1961]  finds  the  mode  sum  for  a  flat  waveguide  and  then 
multiplies  by  a  geometric  factor  to  correct  for  sphericity.  Wait  [1961]  computed  the  attenuation 
and  phase  velocity  of  the  least  attenuated  residue  for  certain  specified  models,  including  in  the 
solution  some  of  the  effects  of  anisotropy,  and  Wait  and  Spies  [1963]  considered  the  amplitude  of 
the  "excitation  factor"  of  the  first  two  residues  with  a  perfectly  conducting  earth  and  ionosphere 
and  of  the  first  residue  for  a  finitely  conducting  ionosphere. 

The  possible  multimode  character  of  the  guide  has  been  noted  by  Wait  [1962a,  b];  however, 
since  large  scale  computer  techniques  have  not  been  available,  httle  has  been  done  to  sum  aU 
modes.  Indeed,  at  distances  less  than  about  3000  km,  it  is  not  sufficient  to  consider  only  a  single 
such  mode,  if  the  frequency  is  greater  than  about  15  kc/s,  for  a  daytime-noon  ionosphere  model 
employed  by  Johler  and  Berry  [1962,  1963]  or  at  frequencies  greater  than  about  4  or  5  kc/s  for  a 
corresponding  nighttime  model  ionosphere.  Johler  and  Berry  [1962,  1963]  computed  LF,  VLF, 
ELF  fields  by  summing  the  classical  series  of  zonal  harmonics  and  by  summing  the  residues  for 


'  This  paper  was  presented  at  the  VLF  Symposiutn,  Boulder,  Colo.,  on  August  12,  1963. 
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the  contour  integral  of  the  Watson  transformation  (called  "modes"  by  Wait  [1962a,  b]  and  by 
others).  The  present  paper  extends  the  previous  work  of  Johler  and  Berry  [1962]  to  take  account 
of  anisotropy.  Curves  are  presented  for  the  field  which  include  all  residues  or  modes  of  signifi- 
cance to  the  computations. 

2.    Theoretical  Solutions 

The  electromagnetic  propagation  of  radio  waves  in  an  anistropic  medium  is  described  by 
Maxwell's  classical  equations.    For  the  continuous  time-harmonic  solution,  with  a  frequency, 

d  ^     .         d  . 

—  E  —  i(oE,  —  H  =  i(oH, 
at  dt 

V  xE  +  fjLo  i(oH  =  0 

V  X  H  -  [e]  •  ^=  0,  (1) 

— >  — > 

where  E  and  H  are  the  electric  and  magnetic  intensities  (rationalized  mks  units)  and  /jlq,  to  are  the 
permeabiUty  and  permittivity  (€o=  /-to  where  c  is  the  speed  of  light)  of  space  and  [e]  is  the  di- 
electric tensor  of  the  anisotropic  medium.    In  the  isotropic,  nonconducting  case. 


[e]  =  eico  =  eoito  in  free  space. 


(2) 


The  dielectric  tensor  in  an  anistropic  medium,  employing  the  Cartesian  coordinate  systems, 
X,  y,  z,  figure  1,  can  be  written  [Johler  and  Harper,  1962;  Johler,  1963a,  1963b] 
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The  motion  of  the  electron  is  described  by  the  Langevin  equation, 

mio^+  mgV+  eiVx  1b)  +  eE  =  0  (4) 


where  the  complex,  frequency-dependent  parameter,  g,  is  assumed  to  be  a  constant,  g  =  v  [Johler, 
1963a]  where^  Wh  =  fJ-o  ej^lm,  the  gyrofrequency  for  electron  charge,  e,  mass,  m,  and  the  magnetic 
field,  J^—  or  the  magnetic  induction  B  =  (Xo^.  0)'^=  Ne'/eom,  the  angular  plasma  frequency 
squared,  /  is  the  magneti^  dip  measured  from  the  horizontal,  A^is  the  electron  density  of  the  plasma 
and  the  magnetic  field,  jCis  assumed  to  be  oriented  in  the  y-z  plane. 
If  a  plane  wave, 

(5) 
where 


E  =  A  exp 


iu)t  ■ 
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T7D  =  x  sin  (/),  sin      +  y  sin  </),  cos  (/)„  +  2^ 
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Figure  1.    Cartesian  coordinate  system  and  spherical 
geometry. 


is  assumed  as  a  solution  of  Maxwell's  equations  in  such  a  medium,  where  (f)i  is  the  angle  of  incidence 
relative  to  the  vertical  z-direction,  and  is  the  magnetic  azimuth  measured  clockwise  from  the 
y-z  plane  containing  the  magnetic  vector,  with  a  dip  angle  /  from  the  horizontal,  a  quartic  in  the 
complex  number,  ^,  is  deduced  [Johler  and  Walters,  1960;  Booker,  1939] 

a4i^^  +  a3^^  +  a2i^'  +  arC'+ao  =  0.  (6) 

The  coefficients  of  this  quartic  are  given  in  appendix  1,  and  a  complex  index  of  refraction,  iq,  of  the 
medium  is  defined, 

Equation  (6)  is  analogous  to  Booker's  [1939]  quartic  in  classical  magneto-ionic  theory,  because 
Booker  first  attached  physical  significance  to  the  four  roots  as:  ordinary  (o)  and  extraordinary  (e), 
upgoing  (i)  and  downgoing  (r)  propagation  components,  17  =  17^,' p.  The  downgoing  components 
can  arise  from  boundaries  (rnore  than  one)  or  electron  density  gradients  in  the  z-direction.  Thus, 
in  the  homogeneous  medium,  only  upgoing  ordinary  and  extraordinary  propagation  components 
exist. 

Attempts  to  solve  Maxwell's  equations  in  an  anisotropic  medium  in  coordinate  systems  other 
than  Cartesian,  such  as  spherical,  r,  6,  (f),  figure  1,  coordinate  systems  lead  to  intractable  equa- 
tions by  conventional  methods.  However,  for  propagation  theory  in  the  space  between  the  plasma 
and  the  terrestrial  sphere,  approximate  solutions  can  be  worked  out  and  tested,  which  task  is  an 
objective  of  this  analysis. 

The  spherical  geometry  of  the  model  propagation  media  is  shown  in  figure  1,  together  with 
the  transmitter,  receiver,  and  boundaries.  The  origin  of  a  spherical  coordinate  system  (r,  6,  (f))  is 
the  center  of  a  homogeneous  terrestrial  sphere  with  a  radius,  a,  and  a  wave  number, 

,      at    /       .  o-A^qC^ 

^.  =  -^e.-.—  (8) 

where  €2  is  the  dielectric  constant  relative  to  a  vacuum,  or  the  permittivity,  6  =  6260  and  cr  is  the 
conductivity  (ohms  "Vmeter)  of  the  terrestrial  sphere.  The  space  between  the  ionosphere  and 
the  terrestrial  sphere  is  represented  electrically  by  the  wave  number, 

^i=-T7.  (9) 
c 


where  tji  is  the  index  of  refraction  (tji  1  for  vacuum,  rji  ~  1.000338  for  air  at  the  surface  of  the 
earth). 


A  vertical  electric  Hertzian  (electrical  point  source)  dipole  is  located  at  9  =  0,  r  =  6(a  <  b  <  d), 
and  the  field  is  to  be  calculated  at  any  point  (r,  6,  cj)),  although  the  specializations  in  this  analysis 
confine   the   computations  to  the  space  {a  <  r  <  d).    Azimuthal  symmetry  will  be  assumed 


2.  1.  Isotropic  Model  Ionosphere 


When  the  earth's  magnetic  field  is  neglected,  the  ionosphere  can  be  characterized  by  the 
isotropic  wave  number 


0}{v  +  ioj) 


(10) 


The  solution  for  this  case  is  well  known  [Watson,  1918] . 

The  vertical  electric  field,  in  particular,  can  be  written,  using  the  Watson  transformation,^ 
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where  P„{cos  B)  is  the  Legendre  function. 
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H[!\^l2\s  the  Hankel  function  of  ith  kind,  order  f  +  |,  and  argument,  z,  and  A +  1/2  is  Bessel's  func- 
tion. i}l\z)  represents  an  outgoing,  radial  elementary  spherical  wave  and  l,[}\z)  represents  an 
ingoing,  radial,  elementary  spherical  wave,  and  ^n{z)  =  \  [^*„^'(2)  +  ^'f  (2)] . 
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where 
and 


2  For  purposes  of  graphical  presentation  in  this  paper,  /o/  =  47r/jUflC  =  3.34  (10"-)  ampere  meters. 
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The  zeros,  Vs,  of  D1-112  =  0,  are  ordered  or  numbered  so  that  |  Im  (vm)  |  <  |  Im  {v„)  |  if  m  <  n,  so 
the  residue  series  is  convergent,  since  Im  (vs)  <  0  and 
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(18) 


provided  Ivsl  >  >  1  and  d  is  not  near  0  or  tt.  Therefore,  Re  {vs)  is  related  to  the  phase  of  a  single 
residue,  5,  and  the  attenuation  rate  is  [Attenuation] «  =  1.364  \  Im  {vs)  \  dB/1000  km. 
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weights  the  significance  of  the  residue. 

The  expression  (11)  is  exact  except  for  a  line  integral 
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which  is  zero  if  only  /?[  1/2  is  an  even  function  of  v  [Berry,  1964a] . 

This  requirement  (R1—112  even)  is  much  less  stringent  than  the  conditions  given  by  Wait  [1963b, 
1962b],  i.e.,  that  the  impedance  boundary  conditions  are  exact  at  the  ionosphere  and  the  earth  or 
that  \k2\  »  ki  and  |A;3|  >>  ^1.  Since  (11)  without  the  integral  becomes  exact  as  the  earth's  con- 
ductivity approaches  infinity,  the  integral  represents  the  field  in  the  medium  a  <  r  <  (/(the  space 
between  the  ionosphere  and  the  terrestrial  sphere)  resulting  from  waves  which  have  escaped  after 
traveling  through  the  medium  r  <  a  (terrestrial  sphere).  The  neglect  of  this  integral  (19)  in  ter- 
restrial radio  wave  propagation  is  therefore  a  very  accurate  approximation,  since  waves  traveling 
in  the  ground  are  very  highly  attenuated  before  again  emerging  into  the  space  a<r<d. 


2.  2.  Anisotropic  Solution  with  Constant  Radial  Magnetic  Field 

A  constant,  radial  magnetic  field  (/  =  =h90°,  (3))  is  now  imposed  on  the  ionosphere  so  that  it  is 
no  longer  isotropic.    The  permittivity  tensor,  (3),  then  becomes, 
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and  where  the  upper  or  lower  sign  (±,  T)  is  taken  for  a  dip  angle  /  — —  90  or  +90°,  respectively 
(the  dip  or  inclination  angle  is  measured  from  the  horizontal  by  convention).  Wait  [1963a]  de- 
rived a  mode  equation  for  this  model  (vertical  magnetic  field,  /  =  +  90°,  i.e.,  directed  in  the  —  f 
direction).  Following  Wait's  procedure  but  writing  the  solution  first  as  a  series  of  zonal  harmonics, 
then  making  the  Watson  transformation,  the  Hertz-Debye  potential  is  found  in  the  space  a<  r<  d. 
The  vertical  electric  field  (11)  now  becomes 
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where  Bv  and  Ev  are  6X6  determinants. 
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Zg  is  the  ground  impedance, 
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in  which  Zo  ~  12077  ohms  is  the  impedance  of  space,  uq  and  Ue  correspond  to  the  ordinary  and 
extraordinary  waves,  either  magnetic  or  electrical  solutions  of  an  assumed  planar  form  in  the 
ionosphere:  A  exp  [— Uo,  p^],  where  A  is  an  appropriate  constant.  By  analogy  with  the  Booker 
[1939]  quartic,  (6),  uq  and  Ue  are  solutions  of 


+  {s\  +  5i)u2  +  s\si  -  si       =  0, 


(28) 


such  that  Re  («o,  e)  >  0  (upgoing  propagation  components),  where 
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si  =  k'i--X\  (29) 
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ao,.  =  —   (32) 

OigUo,e 

(35) 


Equation  (28)  is  comparable  to  (6),  the  Booker  quartic,  when  7  =  90°,  ipa  =  0. 
The  first  two  columns  of  Bv  are, 

-  ClKk,d)C!Mk^b)        i^lKk^r)C^l\k^d)  -  C['Kk^r)Cf{k^d) 

-(,^f-,ikid)l'lKk^b)  C^«(^ir)^<iV'^i^-^';,'(^iOC'!\(A:i(i) 

0  0 
0  0 

0  0  (36) 

and  columns  3  through  6  are  the  same  as  the  corresponding  columns  of  Ec,  (23). 

Equations  (24)  and  (25)  result  from  the  use  of  the  impedance  and  admittance  boundary  condi- 
tions [Wait,  1963a]  ,^ 

He  =  YgE>p 

Ee  =  ~ZgH^       r=a.  (37) 
If  the  usual  boundary  conditions  are  used,      and  Bv  are  8X8  determinants. 

Hr,  He,  and  E^  are  no  longer  zero,  since  a  weak  TE  field  is  excited  by  the  corresponding  TM- 
field  (originally  excited  by  the  vertical  electric  source  current)  at  the  boundary  of  the  anisotropic 
ionosphere.    If  the  magnetic  field  is  zero,  (22)  reduces  to  (11)  if  only 


and 


p,.-x/.--.Vl-(^)  (38) 


k^ 

ia)€oZg  =  ij^^.  (39) 


Equations  (38)  and  (39)  are  good  approximations  for  practical  cases  at  VLF. 


■''  For  more  details  see  Wait  [1963a].  Note  that  the  function  A'^'U)  defined  by  Wait  is  identical  to  $'^'(z).  The  latter  notation  is  used  here  since  it  has  over  50 
years  of  historical  precedence  [March.  1912]  and  has  been  used  by  other  more  recent  authors,  e.g.,  Fok  [  1946]  and  Northover  [1962]. 
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2.  3.  Reflection  Coefficients 


Reflection  coefficients  for  a  plane  ionosphere  with  plane  wave  incidence  have  been  treated 
by  Johler  and  Walters  [1960],  Johler  and  Harper  [1962,  1963],  and  Johler  [1963a].  A  method 
is  now  discussed  in  which  these  Fresnel  type  reflection  coefficients  are  introduced  into  the  spherical 
wave  solutions  of  Maxwell's  equations.  The  x-y  plane  becomes  the  surface  of  the  ionosphere; 
y  corresponds  to  6,  figure  1,  z  to  r  and  x  to  —  0.  The  dielectric  tensor  (3)  is  employed  to  deduce  a 
reflection  process  for  waves  in  the  guide  (a  <  r  <  d). 

Dimensionless  field  ratios  are  defined,  Johler  and  Berry  [1962], 


no  =  ^s—  0'  =^  ^  =  — 


no  ^  _  £j:  ^  £^  ne  , 


/r(e)  r(.e) 


(40) 


(41) 


where  o  and  e  refer  to  ordinary  and  extraordinary,  respectively.  These  ratios  are  functions  of  the 
angle  of  incidence  and  the  connecting  relationship  is 

&m<pi  =  (n  +  y2)lkid.  (42) 
The  fields  which  satisfy  Maxwell's  equations  (1)  in  the  guide  {a  <  r  <  d),  can  now  be  written, 

(a<b<d) 


V  ki  l^,i{kib)Qn'{kir) 

Hr  =  S  [en  &'ikir)  +f„  il,„ik,r)] 

I,   kxjJLQlO)  '- 

E6  =  2.-  w.'v^   +6„dr'('^^i'-)  +  c«0„(/:ir) 

»  ki  lilfuikib)  C>P  (^ir) 

He  =  I.^[e„  i^ik^r)  +/„  ,^,;(A,r)] 


r<b 
r>b 


r<  b 
r>  b 


n  "^1 


B 


n  fJiolCO 


X^Vikrb)  iPnikir) 
iPnik.b)  ^\fHk,r) 


+  bnCn'if^ir)  +  Cnlpn{kir) 


r<b 
r>b, 


(44) 

{a<b<d)  (45) 
(46) 
(47) 

(a<b<d)  (48) 


where 


A  = 


B  = 


n(n  +  l)(2n  +  l)P„(cos  0) 


f'b-' 

sin  6')(2n+l)P'(cos  0) 


hl/J-oC 


rb' 


477- 

hlfJ^oC 


For  the  waves  in  the  ground,  k\  is  replaced  by  k>  (8),  source  terms  containing  6„  are  zero,  and 
in'KkiO)  waves  also  vanish.  The  constant  c«  is  replaced  by  a,,  and  h,„  the  boundary  constants 
for  waves  in  the  ground  for  the  electrical  and  magnetic  solution,  respectively.  Waves  in  the 
ionosphere  are  written,  as  the  product  of  a  field  component  F  and  a  boundary  constant. 


for  ordinary  and  extraordinary  components.  Similarly, 
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H^,n  =  dnF^^}  (49) 

In  the  ionosphere  the  known  solutions  of  the  field  ratios  (40)  and  (41)  are  introduced  to  yield  bound- 
ary conditions  which  satisfy  all  tangential  field  types  which  can  exist  at  the  boundary, 

£"61. 1  =  Ee,  3 

H<p,  I  —  H^,  3 

He. :  =  Qe.nH,.  3  =  ^<  W.'s  +  Qi'rlH!^:,  •  (50) 
These  boundary  conditions  can  be  written  in  matrix  form  using  the  abbreviations      =  [}^\kia). 


^1 


^2 


^1. 


4>ia 

k, 
(//id 

k, 


k, 
^1 


^id 

k, 
A. 


k-z 


a,i 

bn 

Cn 

(—  l/'lft)  l\d 

dn,o 

ki 

d  n,e 

en 

1 

fn 

1 

hn 

In  the  guide  (a<r<d)  the  fields  which  uniquely  satisfy  Maxwell's  equations  (1)  for  the  specifi- 
cation of  the  sources,  can  now  be  determined  with  the  aid  of  the  boundary  constants,  which  can 
be  written 

^U,  +  bn  =  (^W  -  R%U)(l  -  RUlrmWn 

Cn  =  (1//16  -  RlWi-  n,a  -  R%T%,J  -  RmTl^TUDn 
en^{4'ib-RlCib)R',„TlJDn 

fn  =  (llflO  +  Rlill>)njDn 


where, 


Dn      (1      /?|7'|p)(l      R%,T%„^)  R^Rf,^T^„^Tf,^p 


or  in  matrix  notation,  D„  is  determinant  of  the  matrix 


[1] 


7?^,  0 
0  /?? 


•Ts  Ts 
■'■  me      ^  mm  J 


(52) 


(53) 
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The  spherical  type  reflection  coefficients  are, 


C^Kk^a)  k^iSmha) 


+  Rn 


—  1  ^    ^  1  ijj'nikia) 
^  ifjnikia)    ki     "     k-z  ^n(kia) 

^1    "  -^2 


Ts 
^  e 

^  Hi/ 


ipnikid) 

xjjnikid) 
ifjnikid) 


^nikid) 

where,  for  the  case  ipa  =  0  in  plane  wave,  plane  ionosphere  solution. 


(54) 


(55) 


1 

C(Cn) 

0 
0 


-1 


m 

nm 

Vpm 


_0(o) 


-1 


me 


0 
0 

1 


(56) 


1 

Ciiljjn) 

0 
0 


0(0) 
I   •  /Vme 

'  Vem 

-1 


I    •  /Vme 
'  Vem 

-1 


-Q(e) 
Vme 


1 

Ci(C«) 

0 
0 


(57) 


Tmm  ^ 


1 

C,(i//„) 

0 
0 


I    •  /Vme 
'  Vem 

-1 


)(0) 

'me 


)(e) 


Ve 


irjeV^I^' 


)(e) 

me 


0 
0 

C,(^«) 

1 


(58) 


and, 


rme  =  A- 


A  = 


0 
0 

CiiCn) 
1 

1 

Ciiljjn) 

0 
0 


I    •  /V  me 

'  Vem 

-1 
Vme 


Q(0) 
I    •  /Vme 

'  Vem 

-1 


)(0) 
'me 


+  iT?e 


(e) 
em 


Vme 


+  tT7e 


)(e) 
'me 

n(e) 
Ve»n 


)(e) 
'me 


0 
0 

Ciiljjn) 
1 

0 
0 

C#n) 
1 


(59) 
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Ciil^n)  =  l,  ~  COS  4>i 


cos  4>i 


{ktd}^ 


(60) 
(61) 


and  in  the  Watson  transformation,  n  =  fs-1/2. 

A  vertical  magnetic  dipole  or  an  elementary  line  magnetic  current,  Iml,  volts-meters,  is  equiva- 
lent to  an  elementary  current  loop  loO,  ampere-(meters)^,  w^here  O  is  the  area  of  the  loop  in  square 
meters. 

Iml^  /J-oifoIoO. 

The  magnetic  dipole  source  solution  can  be  obtained  from  (43)  to  (48)  by  moving  the  source  func- 
tions ^ib^tr,  ^ib^i'r  ,  or  ^ir'/'ift,  ^I'r  i/'i6  from  (43),  (45),  and  (48)  to  (44),  (46),  and  (47)  respectively,  and 
replacing  lol  by  ikJoO.  The  fields  propagated  by  such  a  vertical  magnetic  dipole  are  also  equiva- 
lent to  the  fields  propagated  in  the  equatorial  plane  of  a  horizontal  electric  dipole. 

The  Watson  transformation  can  be  introduced  to  facilitate  computation,  which  for  the  par- 
ticular case  of  the  vertical  electric,  Er,  field  can  be  written. 


4-77 


lol 


-Itt 


r[  cos  {vsTt) 


\'s-l/2(— cos  0)  X 


Ks-idl^s-U-iikir)  +  c[,^-ii2^vs-il2(kir) 


(62) 


dv 


D1.-1I2 


where  6'.  =  Dibv,  c'v  =  CiDv,  and  the  poles  v  —  Vg  are  the  roots  of  the  mode  equation  (53), 

Dvs-ll2  =  0. 

Computations  can  be  made  more  convenient  by  introducing  the  exact  relationship 

Wz)  =  He'W  +  C'(2)] 


(63) 
(64) 


into  (51).  The  result  is  that  ijjnix)  functions  are  replaced  by  ^*^'(^)  functions  in  D„  and  in  the  defi- 
nitions of      and      (54),  (55),  and  the  quantity  in  square  brackets  in  (43)  is  (assuming  r  =  a  =  b) 


y(2) 


—  ffs  rwri^)  in  —Ts   j?s  vi         Ts  \ j's  Ts  jfs  n-1 

>-(l)    e  /  felaSla  I  V-^      ^  mm^^mJ\^      y{2)    ee>      >.(2) -*  em-*  me"-m  /  ■'^n 


which  agrees  with  the  result  obtained  by  Wait  [1963b]  using  impedance  boundary  conditions 
instead  of  (50).    Thus,  the  impedance  boundary  conditions  [Wait,  1963b]  are. 


'Ee 

Zu 

-Zri 

He 

.-£'<*. 

_  Z21 

Z22. 

Z\2 

E,_ 

Z2^ 

Z22  _ 

where  the  impedance  Zy  must  be  expressed  with  reflection  coefficients  and  the  cosine  of  the  angle 
of  incidence  Ct, 

Zn  =  —  2  TmeZolD 

Zi2=  — CiZo[(l  —  Tmm)ii  ~  Tee)~  TemTme]ID 
Z21  =Zo[(l  +  7'mm)(l  +  Tee)  —  TemTmeVDCi 
Z22—     2  TemZ o/D 

where  D  —  {1  —  Tmm){i  +  Tee)  +  TemTme-  It  is,  therefore,  concluded  that  the  impedance  boundary 
conditions  and  the  conditions  stated  (50)  are  equivalent. 

The  analysis  leading  to  (43)  to  (48)  and  (62)  is  clear  for  propagation  in  the  guide  [a  <  r  <  d) 

subject  to  some  restrictions,  such  as  ^  =  0,  4>a  =  0,  180°,  /  — ±90°.    For  cases       7^  0,  (ha^O, 

o(p  d(p 
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180°,  I  ¥=±  90°,  the  approximation  accuracy  can  be  tested  by  a  summation  of  the  spherical  har- 
monics [Wait,  1963b]  in  the  variable  </>. 

3.    Computational  Methods 

The  solutions  in  section  2  were  programmed  for  a  large-scale  digital  computer.  This  allows 
accurate  computation  of  the  functions  involved,  makes  summation  of  all  significant  modes  feasible, 
and  provides  for  study  of  the  field  variation  with  any  parameter. 

The  functions  are  computed  using  (12)  and  [Fok,  1946] 


where  ^^^(tanh  w  —  w),  and  cosh  w  =  vlz.  The  approximation  for  H[!'Hz)  involving  the  Airy  inte- 
gral function,  Wk{z),  which  was  introduced  by  Fok  [1946],  is  the  first  term  of  a  series  in  ^1— 

representing  the  right  side  of  (65),  so  it  is  equivalent  to  (65)  when  v  =  z.  However,  (65)  is  uniform 
over  a  wide  range  of  v  and  z.  Berry  [1964b]  gives  detailed  instructions  for  computing  (65)  and 
finds  the  magnitude  of  the  relative  error  is  about  10"^/|z|. 

Taking  the  derivative  with  respect  to  order,  one  finds  that 


A        .  ^  exp  [(-  1)^-+ Vi/6] 


+  iw  exp  [(-  l)''+'27Til3]m%  (-  iO 


(66) 


The  determinates  are  computed  numerically  by  reducing  them  to  equivalent  triangular  de- 
terminates using  the  Crout  forward  solution  and  then  taking  the  product  of  the  diagonal  elements 
[Hildebrand,  1946].    Full  advantage  is  taken  of  the  zeroes  in  the  corners. 

Let  Dniv)  be  an  n  X  n  determinate  whose  elements  are  functions  of  v,  and  let  D^'\v)  be  the  n  X  n 
determinate  which  is  identical  to  Dniv)  except  that  the  elements  of  the  tth  row  (or  column)  are 
replaced  by  their  derivatives  with  respect  to  v.  Then 

^  Dniv)  =  f^D^i^  ivy  (67) 

^  i=l 

Equation  (67)  is  used  to  evaluate  the  denominator  of  the  excitation  factor,  and  in  (68)  below. 
Given  a  A:th  approximation,  1;^',  to  a  zero,  Vs  of  Dv-nt,  a  (A:+  1)  approximation  is  obtained  with 
Newton's  iteration  [Hildebrand,  1956] , 


+  1)  =:  „(/c)  !^J/2_ 

dv  "s-"^ 


The  iteration  converges  quadraticaUy  to  Vs  if  vi'^'Hs  sufficiently  close  to  Vs. 

Initial  approximations  to  the  Vs  are  obtained  several  ways.  Contours  of  \Dv-ii2\  can  be  drawn 
in  the  complex  f-plane  and  the  zeros  approximated  graphically.  This  is  the  most  foolproof  way  of 
insuring  that  all  significant  zeros  are  found,  but  is  slow,  tedious,  and  expensive. 

For  frequencies  less  than  about  15  kc/s,  kia  is  an  adequate  first  guess  for  fi,  and  first  guesses 
for  higher-order  modes  for  the  isotropic  case  can  be  obtained  from  analytic  expressions  given  by 
Wait  [1957]. 

For  frequencies  above  15  kc/s,  a  first  guess  of  kia  frequently  converges  to  v-z- 
To  find  the  zeros  in  the  anisotropic  case  the  zeros  for  the  corresponding  isotropic  case  are 
found  for  both  vertically  and  horizontally  polarized  source  dipoles  (TM  and  TE  modes).    These  are 
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usually  adequate  first  approximations  for  the  anisotropic  case,  as  shown  in  table  1,  which  gives  the 
zeros  for  /=  10  kc/s,  h  =  70km,  N  =  630,  and  v  =  W. 


Table  1 


•  Isotropic 
TM  modes 

Anisotropic 
modes 

Isotropic 
TE  modes 

1  331   K.Ci      a  (\A  .' 
lOol.oy  —  o.U'+  I 

1313.68-2.14  / 

1313.17-1.79  i' 

1260.09-22.49  / 

1262.34-21.4  i 
1222.85-8.64  i 

1220.62-7.07  / 

1079.42-48.77  / 

1077.35-49.5  i' 
10.56.52-23.9  i 

1048. .54-  18.43  i 

760.67-70.60  / 

739.34-83.7  i 

If  the  Vs  are  known  for  some  set  of  parameters,  they  can  be  used  as  first  approximations  for  a 
case  in  which  one  (or  more)  of  the  parameters  has  been  changed  by  a  small  increment.  The  param- 
eter(s)  can  be  incremented  again  and  again  until  the  desired  set  is  reached.  When  this  is  done 
with  only  one  parameter,  the  result  is  the  variation  of  the  field  with  that  parameter,  which  is  often 
of  interest.    (See,  for  example,  figs.  2,  3,  and  11.) 

The  Legendre  function,  Pv(—  cos  d}  =  Pv[cos  (tt  — 0)]  is  computed  from  [Fok,  1946] 


Pv(cos  6)  = 


r(t;+l) 


Y{v  +  3/2)  V27r  sin  d 


[  e'-^G, + :/2(  1  -  x)  +  e  -       + i/2(x)] , 


(69) 


irhere 


and 


4)  =  e{,v  +  h)-TTi^,x  = 


le" 


2  sin  0  ' 


■x  + 


p.32 


221!(i;+l)       242!(2;+l)(i;  +  2) 


+ 


p.32 


(2n-l)2 


2^"nliv+l)(v  +  2)  .  .  .  {v  +  n 


-x"  + 


(70) 


The  series  dix)  converges  if  v  is  not  a  negative  integer  and  |:x:|<l,  i.e.,  3O°<0<15O°.  If  6  is  outside 
this  interval,  (69)  is  an  asymptotic  representation.  If  only  the  first  term  of  (70)  is  used,  (69)  reduces 
to  (18),  and  the  relative  error  is  one  percent  or  less  if  \v  sin  d\  >  10.  Equation  (69)  is  adequate  at 
VLF  for  all  points  more  than  50  km  from  the  source  and  antipode  and  can  be  used  down  to  15  km 
at  30  kc/s. 

When  6  is  near  tt  [Magnus  and  Oberhettinger,  1949], 


P.(-cos  e)  =  Jo(cl>)  +  sm'  + 6 


■Acf>) 


O 


sin'' 


2 


(71) 


where  0  =  (2v  +  1)  sin 
source. 


TT-e 


Equation  (71)  is  useful  for  studying  fields  near  the  antipode  of  the 

4.  Discussion 


The  computer  program  described  in  section  3  was  used  to  compute  the  field  for  a  variety  of 
parameters.    Three  models  of  the  lower  ionosphere  (sharply  bounded  type)  were  used  extensively: 

Model  I:  h  =  67.5  km,  N  =  56  electrons/cm^,  v=1.6{  W); 
Model  II:  /i  =  75km,/V=160electrons/cm3,  i/  =  4.5(10«);  and 
Model  III:  A  =  70  km,  A^  =  630  electrons/cm^,  10^. 

Model  III  is  a  model  used  by  Wait  [1957]  (corresponding  to  (o,  =  2(10'')  in  Wait  paper).  In  his  early 
work,  Wait  [1957]  used  the  Debye  approximation  [Debye,  1910;  March,  1912] 
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for  all  the  values  of  v  and  z.  Equation  (72)  is  inadequate  at  frequencies  above  approximately  10  kc/s 
for  the  first  or  most  important  mode  (using  Wait's  terminology).    It  is  interesting  to  note  that  the 

inadequacy  of  this  approximation  was  discussed  by  March  [1912]  for  sin  c^,  =— ^ —  =  1,  yet  he  used 

the  approximation.  Also,  von  Rybczynski  [1913]  noted  the  more  accurate  Hankel  approxima- 
tion (65)  and  used  this  as  a  basis  for  more  accurate  corrections  to  March's  [1912]  wrork.  However, 
for  z  large  and  v  small  or  v  large  and  z  small  the  approximation  is  quite  excellent.  The  use  of 
the  approximation  (72)  is  imphed  by  the  use  of  plane  reflection  coefficients  at  the  ionosphere 

m)ii-  ^  ^'  ^^^^  However,  in  this  case  it  is  a  comparatively  simple  matter  to 

replace  cos  0;  with  the"  functions,  Cn'i^id)ICn{kid),  instead  of  the  approximation  real  cos  or 


complex  cos  0;  =  ,  / 1  _  /  ' 
V       \  kid  J 


More  recently,  Wait  [1961]  and  Spies  and  Wait  [1961]  used  an  approximation  containing 
Airy  integrals  which  is  nearly  equivalent  to  (65).  The  characteristics  of  the  least  attenuated  mode 
were  calculated.  This  is  quite  accurate  [Berry,  1964a].  However,  one  mode  is  not  adequate  to  de- 
scribe the  complete  field  which  satisfies  Maxwell's  equations  at  either  intermediate  or  great  dis- 
tances and  higher  frequencies,  figures  2,  4,  5,  and  6,  for  daytime  models  of  the  lower  ionosphere 
(highly  absorbing  models  such  as  model  I).  In  the  case  of  a  nighttime  model  even  more  modes 
are  needed  at  all  distances,  figure  2.  Since  the  prime  objective  of  the  analysis  of  this  paper  is  to 
find  the  field  indicated  by  Maxwell's  equations,  and  develop  accurate  techniques  to  calculate  such 
a  field  rather  than  attempt  to  study  the  physical  significance  of  modes,  all  modes  of  significance 
to  the  calculation  were  summed  in  an  ordered  series. 

The  attenuation  rate  of  the  more  important  modes  as  a  function  of  the  electron  density,  A'^,  in 
the  homogeneous,  isotropic  model  of  the  ionosphere  is  shown,  figure  2,  for  25,  10,  4,  and  1  kc/s. 
Notice  that,  for  a  fixed  frequency,  a  given  attenuation  rate  can  be  found  for  as  many  as  three  dif- 
ferent values  of  N,  whose  magnitudes  are  very  roughly  in  the  proportion  1:10:10^.  This  explains 
why  early  investigators  [Watson,  1918,  1919;  Bremmer,  1949]  deduced  a  high  conductivity  for 
the  ionosphere,  and  why  values  of  A'^  as  different  as  those  in  models  I  and  III  have  been  considered 
more  recently. 

At  VLF,  an  increase  in  electron  density  can  cause  an  increase  in  attenuation  rate  for  model  I 
(daytime-noon  model)  while  at  night  the  attenuation  rate  is  decreased,  figure  2.  The  broad,  nearly 
flat  minimum  in  the  attenuation  curves,  together  with  the  insensitivity  to  a  magnetic  field  pre- 
viously mentioned,  helps  explain  the  observed  stable  characteristics  of  VLF  propagation  in  the 
vicinity  of  10  kc/s  and  indicates  that  accurate  determination  of  ionospheric  characteristics  using 
VLF  radio  waves  may  require  sophisticated  techniques.  The  attenuation  rate  seems  to  be  most 
sensitive  to  the  assumed  height  of  the  ionosphere,  figure  3.  This  insensitivity  to  the  detailed  char- 
acteristics of  the  lower  ionosphere  suggests  use  of  higher  frequencies,  such  as  LF  (100  kc/s)  which 
are  more  sensitive  to  ionosphere  changes  rather  than  VLF  (10  kc/s),  to  study  experimentally  and 
theoretically  the  detailed  characteristics  of  the  lower  ionosphere. 

Figures  4,  5,  and  6  illustrate  the  amplitude  of  the  field  and  the  phase  correction  or  phase  lag, 
(pc=|  ^TgEr  \  —kiad,  of  the  vertical  electric  field  as  a  function  of  distance  from  the  source,  for 
30,  20,  6,  1,  and  0.1  kc/s,  for  three  models  of  the  ionosphere.  At  30  kc/s  the  field  shows  large 
oscillations,  even  at  7000  km.  The  attenuation  of  the  second  mode  is  nearly  as  small  as  that  of 
the  first  mode  at  this  frequency,  while  the  excitation  factor  of  the  first  mode  is  Vs  to  Vio  that  of 
the  second  mode.  Thus,  the  two  modes  interfere  with  each  other,  even  at  great  distances. 
The  decrease  in  the  excitation  factor  of  the  first  mode  (s—  1)  with  increasing  frequency  was  noted 
by  Wait  [1962b]. 

At  the  point  at  which  the  curve  for  the  amplitude  of  the  field  as  a  function  of  distance  shows 
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Figure  3.  Attenuation  rate  of  first  three  modes  as  a 
function  of  altitude,  h,  of  model  ionosphere  boundary 
at  10  kcjs. 
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Figure  4.    Amplitude  and  phase  of  the  vertical  electric 
field,  Er,  as  a  function  of  distance  for  model  I. 
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Figure  6.    Amplitude  and  phase  of  the  vertical  electric 
field,  Er,  as  a  function  of  distance  for  model  111. 
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Figure  5.    Amplitude  and  phase  of  the  vertical  electric 
field,  Ef,  as  a  function  of  distance  for  model  II. 
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Figure  7.    Amplitude  and  phase  correction  of  the  vertical 
electric  field  as  a  function  of  frequency  for  model  I. 


Figure  8.    Amplitude  and  phase  correction  of  the  vertical  FIGURE  9.    Amplitude  and  phase  correction  of  the  vertical 

electric  field  as  a  function  of  frequency  for  model  II.  electric  field  as  a  function  of frequency  for  model  III. 

a  sharp  cusp,  there  is  an  abrupt  change  of  approximately  tt  radians  in  the  phase  correction,  (f)c. 
This  occurs  in  the  20  kc/s  case,  figure  4,  at  approximately  800  km,  and  in  the  30  kc/s  case,  figure  6, 
at  approximately  5300  km.  The  cusps  occur  when  the  first  two  modes  are  exactly  out  of  phase 
and  have  almost  the  same  amplitude. 

Figures  7,  8,  and  9  show  the  amplitude  of  the  field  and  the  phase  correction  as  a  function  of 
frequency  at  2000  and  4000  km  for  the  three  models  of  the  lower  ionosphere  both  with  and  without 
a  vertical  magnetic  field.  All  three  models  show  an  absorption  band,  the  position  of  which  can 
be  changed  by  altering  the  assumed  electron  density,  for  example.  The  presence  of  the  vertical 
magnetic  field  has  a  small  effect  in  model  I,  figure  7,  and  very  little  effect  at  VLF  in  all  models. 
However,  figures  8  and  9  show  a  decrease  in  the  field  at  ELF  when  the  vertical  magnetic  field  is 
present.  Wait  and  Carter  [1960]  noted  this  effect  and  stated  that  it  would  occur  when  the  col- 
Hsion  frequency,  v,  is  comparable  with  the  gyrofrequency  of  the  electron  gas.  Model  II,  figure  8, 
most  nearly  satisfies  such  a  requirement  and  indeed  shows  the  most  decrease  in  the  amplitude  of 
the  field.  The  presence  of  the  static  magnetic  field  tends  to  increase  the  amplitude  of  the  field 
at  VLF,  figures  7,  8,  and  9.  Also,  in  figures  8  and  9,  the  magnetic  field  causes  a  "step"  in  the 
amphtude  curve  as  the  amplitude  of  the  field  decreases  into  the  absorption  region  from  higher 
frequencies,  and  the  minimum  field  strength  occurs  at  a  lower  frequency.  The  "step"  in  the 
curve  is  caused  by  interference  between  the  first  TE-mode  (transverse  electric)  and  the  first  TM- 
mode  (transverse  magnetic).  Ordinarily,  the  excitation  of  the  first  TE-mode  is  approximately 
0.01  that  of  the  first  TM-mode.  Physically,  this  is  understandable,  because  the  TM-  or  normal 
modes  are  excited  by  the  vertical  dipole  sources,  while  the  TE-  or  abnormal  modes  are  coupled 
in  at  the  anisotropic  ionosphere  boundary.  However,  at  about  3  or  4  kc/s  in  models  II  and  III, 
figures  8  and  9,  the  excitation  of  the  first  TE-mode  was  comparable  with  that  of  the  first  TM-mode, 


17 


and  its  attenuation  rate  was  greater,  so  that  they  interfered  with  each  other.  Figure  10  illustrates 
the  amphtude  of  the  field  as  a  function  of  distance  for  two  such  cases.  The  change  in  the  slope 
of  the  curve  clearly  indicates  the  change  in  normal  or  abnormal  mode  dominance. 

The  effect  of  arbitrary  dip  or  incUnation,  /,  of  the  magnetic  held  vector,  on  amplitude  of  the 
field  is  shown,  figure  11,  as  a  furlction  of  dip  angle  at  10,  1,  and  0.1  kc/s.  The  change  in  the  ampli- 
tude of  the  field  with  dip  angle  is  more  obvious  at  1  kc/s.  Comparison  of  figures  6  and  11  shows 
dip  angles  between  30  and  60  deg  exhibit  virtually  no  difference  between  the  isotropic  and  aniso- 
tropic curves  at  10  kc/s  (VLF). 

As  a  consequence  of  the  numerical  and  analytical  insight  gained  by  the  above  described  analy- 
sis, the  procedure  to  foUow  in  introducing  reflection  coefficients  is  evident.  Considerable  caution 
is  necessary  since  the  Debye  type  approximation  in  the  guide  {a  <  r  <  d)  (64)  is  again  not  accurate 
for  reasons  described  above.  But  in  this  case,  there  is  httle  difficulty  calculating  the  spherical 
wave  functions  Cn^{z)  and  Cn^'iz). 

The  analysis  presented  in  this  paper  has  been  a  theoretical-mathematical  study  based  on  physi- 
cal premises  and  a  number  of  phenomena  have  been  predicted.  Verification  of  these  phenomena 
have  not  as  yet  been  made  under  carefully  controlled  experiment. 

5.  Conclusions 

The  effect  of  anisotropy  in  the  ionosphere  on  the  field  of  LF,  VLF,  ELF  radio  waves  can  be 
described  theoretically  for  the  region  between  the  ionosphere  and  the  earth  by  following  conven- 
tional isotropic  analysis  procedures  and  introducing  the  magneto-ionic  theory  with  the  aid  of 
reflection  coefficients.  In  addition  to  a  considerable  number  of  normal  modes,  abnormal  modes 
must  also  be  summed  to  completely  describe  the  field.  The  LF  portion  of  the  long  wavelength 
spectrum  is  more  sensitive  to  ionosphere  modes  and  changes  in  electron  density,  magnetic  field, 
etc.,  than  the  VLF.  The  most  important  parameter  at  VLF  is  apparently  the  height  of  the  model 
ionosphere  boundary  for  the  cases  studied  in  this  paper.  The  results  of  this  paper  suggests  the 
work  be  extended  by  introducing  more  accurate  reflection  coefficients  which  are  based  on  an 
ionosphere  in  which  the  electron  density  varies  with  altitude  and  distance. 


The  authors  are  indebted  to  their  colleague,  James  R.  Wait,  for  his  critical  review  of  this  work. 
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The  solution  of  Maxwell's  equation  (1)  for  plane  waves  in  a  Cartesian  coordinate  system  con- 
sidered locally  at  the  model  ionosphere  [Johler  and  Berry,  1962]  and  [Johler  and  Harper,  1962] , 
for  the  dimensionless  field  ratios. 
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in  which  it  has  been  assigned  that  the  solution  of  equations  (1)  and  (4)  in  Cartesian  coordinates 
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respectively,  where  the  simultaneous  solution  of  equations  (1)  and  (4)  in  Cartesian  coordinates 
yields  upon  ehminating  V  and  E, 
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All  of  the  equations  and  text  of  this  publication  were  composed 
on  the  keyboard  of  a  photocomposition  system,  eliminating  all 
hand  insertion  of  equation  characters. 


